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This paper considers the application of a digital computer and discrete 
Fourier transform (DFT) techniques to the measurement of signals known 
to comprise only single-frequency tones. We discuss the use of weighting 
functions to improve the effective selectivity of a measurement system that 
estimates the frequencies and levels of tones from the coefficients of their 
DFT. We present three classes of weighting functions which may be used 
to improve the inherent accuracy of such a system. The form of the weighting 
functions was chosen to minimize the amount of computer memory required, 
without using too much computer time. Several formulas are derived for 
estimating the frequency and level of a tone from its DFT coefficients. We 
chose the formulas to minimize computation time. 

Simulation results indicate that, through the use of a proper weighting 
function, a DFT measurement system that uses 512 samples taken at a 
sampling frequency of 7040 Hz can be designed so that the maximum error 
in the frequency estimates of two tones near 1000 Hz and separated by 
approximately 50 Hz is about 0.03 Hz. The corresponding maximum error 
in the level estimate is on the order of 0.03 dB. 

I. INTRODUCTION 

There have been numerous articles, in recent years, dealing with the 
use of the discrete Fourier transform (DFT) in the area of spectrum 
analysis. Much of this interest was motivated by the availability of a 
computational algorithm that facilitates the rapid computation of 
DFT coefficients by a digital computer. The algorithm is, of course, the 
fast Fourier transform (FFT). 

We are concerned with the problem of applying DFT techniques to 
the measurement of the levels and frequencies of single-frequency tones, 
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Fig. 1 — A DFT measurement system. 

particularly tones from a data set during a test. Figure 1 shows the 
system we have in mind. A band-limited received signal, known to 
comprise one or more single-frequency tones, is periodically sampled by 
an A-D converter. A total of N samples are taken and the DFT co- 
efficients are computed from the samples. The computer determines 
which of the DFT coefficients are "large", indicating the approximate 
frequencies of the received tones, and then proceeds to compute accurate 
estimates of the frequencies and levels. Methods for achieving the first 
part of the procedure are well known. This paper is devoted to a con- 
sideration of how best to go about the last step in the process, the ac- 
curate estimation of the frequencies and levels of the received tones. 

In data set testing, the tone measurement system would be used 
occasionally during a test and would have to consume a minimum 
amount of real time. Thus we have directed our attention toward estima- 
tion methods that use simple formulas and require a minimum amount 
of computer memory. 

Our attention is confined to the problem of leakage, its reduction by 
smoothing (windowing) functions, and the development of formulas 
which extract tone levels and frequencies from the list of DFT coeffi- 
cients. We don't discuss the important, but secondary, problems of 
round-off errors and other noise sources. 

II. REVIEW OF DISCRETE FOURIER TRANSFORM 

The definition and properties of the discrete Fourier transform are 
discussed in Refs. 1 and 2. The following review is to refresh the reader's 
memory and establish the notation that we will use later. 

2.1 Definition of Discrete Fourier Transform 

Consider an ordered set of numbers {X„} where n = 0, 1, 2, ■ ■ • , 
N — 1. Following Cochran, and others, 1 we define the discrete Fourier 
transform (DFT) of the set {X n } to be another set of numbers, {A K }, 
with 



A K = £ x„ 



-i2mK/N 



all integer K. 



(1) 
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The inverse transformation is 



X n = jj E Agf"*", n - 0, 1, 2, • • • , N - 1. (2) 

2.2 Useful Properties 

Several properties of the DFT are utilized in later parts of this paper. 
The important properties are recorded in this section for future ref- 
erence. Reference 2 provides a more complete list. Derivations are 
included only for results that may not be well known. 

From equation (1) it is obvious that if the X n are real, then 

A_ K = A* K (* denotes conjugate), (3) 

A K+N = A K , (4) 



and 



A„. K = A. K = A* . (5) 



2.2.1 Convolution 
Let 



B K = £ X n e-'*' nK/N (6) 



and 



then 



N-l 



C K = Z Y n e-»* nK/ \ (7) 



A m = E X n Y n e~^ nm/N = ^ E B K C m . K . (8) 

„-0 iV K-0 

In other words, if \B K \ and {C K \ are the DFT of {X n } and \Y n \, 
respectively, then the DFT of {X n Y n } is given by equation (8). 

2.2.2 Power 

It can easily be shown, for X n and A K denned by equations (1) and 
(2), that 

, N-\ N-l 

± E A*Al = E XI . (9) 

iV K-0 n=0 

If the X„ are samples of some function, f(t); that is, if X\ = f(nT/N), 
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then 

rp If -I r-T 

ff-too *» n = J 

if the integral exists. Thus, for large iV, 

f md*«|x;^. do) 

Hence, from equation (9), 

Tf'm*~&%\**r- od 

2.3 Relationship to Fourier Transform 

The DFT of samples of a signal has a simple relationship to the 
regular Fourier transform of the signal. It is instructive to examine this 
relationship. 1 " 

Let g(f) be an arbitrary function, zero for t < and t > T and con- 
tinuous over < t < T. The function is allowed to be discontinuous 
at t = and at t = T. Assume that g(Q+) and g{T—) exist. 

A well-known application of the Poisson sum formula gives 4 

MW + MT-) + %,!§) -%±G<?f) (12) 

where 

G(u) = f g(t)e- iul dt. (13) 

•'0 

Adopting a notation similar to that of Papoulis, 4 we define the "#" 
operation by 

G'(u) = £ E <?(« - &».), ( 14 ) 

i K=-oo 

where 

w. = 27riV/T. (15) 

Then equation (12) can be rearranged to give 

E g($f) = G'(0) + Mff(0+) - (7(7'-)], (16) 

where #(0) is taken to equal #(0+). 

t The recent article by Bergland touches upon this subject and also contains 
an extensive list of references. 3 
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Let h(t) be any function of the sort used above for g(t) with the ad- 
ditional property: h(0) = h(0+). 

Let s(t) be the signal to be analyzed and define 

/(*) = s(t)h(t). (17) 

Let g(t) = i(t)e~' ut and define 

A(o>)=Zj(jf)e- inuT/N . (18) 

Then from equation (16) and the definition in equation (14) we have 

A(») = F'(u) + M/(0+) - KT-)e' iuT ], (19) 

where 

F( u ) = f j{t)e~ iat dt. (20) 

•>0 

If X n = ftnT/N) then the A K defined by equation (1) are given by 

A, = A^f). (21) 

Thus the DFT of the set \j(nT/N) } are points along the curve described 
by equation (19). These points are \/T Hz apart. 

Observe that at w = 2wK/T the term in brackets in equation (19) 
becomes J[/(0+) — f(T— )] which is independent of K and vanishes 
if/(0+) =/(T-). 

2.4 Weighting Functions 

If the DFT is to be taken of the set {s(nT/N)\ for n = through 
N — 1, then h(t) must be a function whose value is unity at t = nT/N; 
n = 0, 1, • • • , N — 1. The function with this property that is usually 
taken to be h(t) is the function h T (t); 

hM _ Jl. 0^t<T; (22) 



[0 , otherwise . 

Other weighting functions, h(t), are often formed by multiplying h T (t) 
by a nontime-limited function. Weighting functions play a very im- 
portant role in systems that use the DFT. The following paragraphs 
attempt to develop and present some of the pertinent theory. 

From equation (19) we see the role that F*(co) plays in A(u>). Since 
f(t) = s(i)h(t), 

F(u) = S(u)*H(u), (23) 
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where the * denotes convolution and <S(a>) and H(co) are the Fourier 
transforms of s(t) and h(t). It can be shown that, subject to the usual 
convergence constraints, 

F'(a) = [S(a)*H(u)]' = S(a>)*H"(a>). (24) 

Thus H(a)), or equivalently H*(u>), plays a central role in the DFT of 
(weighted) samples of s(t). From the development that led to equations 
(19) and (21), we see that, if h(0+) = h(T-), the DFT of samples 
of h(t) is a set of points taken along the periodic curve described by 
H*((j>). It follows, therefore, that the values of h(nT/N) can be obtained 
from 

if) = },£«■&>- 

Also, 

H*(u>) = Z hlfij-y*"*'" - **(0+)[l - •-'•*]. (26) 

Weighting in the time domain is actually done at the points t = 
nT/N; n — 0, 1, 2, • • • , N — 1. For every set of weights to be applied 
at these points there exists a continuous function with the same values 
at the indicated time points. Thus there is no loss of generality due to 
discussing weighting in terms of weighting functions, h(t), that are 
continuous over (0, T) and zero outside that interval. We have to re- 
member, however, that if the set \h(nT/N)} is specified, h(t) is not 
unique. Thus, if H*(u>) is given, h{nT/N) is given by equation (25), but 
h(t) and H(a)) are not uniquely defined. 

There is apparently some confusion in the literature about whether 
H(oi) or H"(o)) is called a weighting function (or windowing function). 
Blackman and Tukey, 5 for example, discuss h(t) and H(w), but when 
Helms 6 writes about weighting with a Dolph-Chebychev function, he 
is evidently referring to H*(u>). More will be said about this later. 
Bingham, and others, in writing about data windows (See Reference 7, 
Part VII) mean h(t). 

Observe that H*(w) is always periodic with period w, , while H{w) 
is not periodic. (If it were, H"(ui) would not converge properly.) Gen- 
erally the H*(u) that one uses will have a prominent main lobe about 
a) = Km, (K is any integer, including zero) and many side lobes. For 
our purposes it is important to obtain a narrow main lobe and low- 
amplitude side lobes. 

The class of H"{u>) with the minimum main-lobe width for a given 
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side-lobe amplitude is known as the (discrete) Dolph-Chebyshev weight- 
ing functions. 8 A convenient form, similar to the one given by Helms, 
but changed to describe the result of weighting by sample values that 
peak at T/2 and are adjusted to cover approximately unit area as later 
weighting functions will do, is the following: 

H"(u>) = | e~ iuT/2 cos [^ cos" 1 [z o cos ^)] (27) 

where the side-lobe amplitude, 1/R, is related to Z by 

R = cosh (N cosh-%) (28) 

and N is the same as used in equation (1). 

The class of H(u) with the minimum main-lobe width for a given 
side-lobe amplitude is known as the continuous Dolph-Chebyshev 
functions, 9 which are unrealizable. The Taylor approximations to the 
continuous Dolph-Chebyshev functions 9-11 are realizable, however, and 
provide almost the same main lobe width for a given maximum side-lobe 
amplitude. 

The problem of choosing "good" shapes for H*(oi) can be approached 
by treating H* (co) or by treating H(w) . Most of the well-known weighting 
functions are discussed in terms of H(ui) or h(t). 

2.5 A Generalization 

If h(t) is a function that is zero for t < T, and t > T, then it can be 
shown (sampling theorem) that i/(co) is given by 

H(u) = re-'"™ sin (0,772) £ t ° n (29) 

„ = _«, 0)1 



2 ~ n7r 



and 

TC n = 



n = H^f). (30) 



Thus the specification of a weighting function is equivalent to the 
specification of the constants, C„ . 



III. SELECTED WEIGHTING FUNCTIONS 



3.1 Leakage and Aliasing 

Leakage will be used here to refer to the problem of the values of 
A (co) due to cos (u Q t + 6 ) interf erring with the values of A (co) at some 
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other frequency, say o>, , where the response due to cos (o>,£ + 0,) is 
to be examined. Leakage, in our system, is minimized by the use of 
weighting functions. 

Aliasing refers to the fact that in a sampled-data system tones with 
frequencies above w,/2 cannot be distinguished from tones with fre- 
quencies less than co,/2. In our system aliasing is avoided by the use 
of the low-pass filter (Fig. 1). 

3.2 Convolution of Weighting Functions 

The object of weighting is to produce the DFT of a weighted set of 
samples of the signal undergoing measurement, s(t). Thus we seek to 
compute 

A K = J2 s(ntMnQe~ i2 ' nK/N , for all K; (31) 

n = 

where t, = T/N. A convenient way of doing the weighting is to first 
compute 

B K =Zs(nt a )e-^ nK/N , (32) 

n=0 

for £ K £ N - 1. Then if the set \H m \, 

II m = E *(n«.)e-"— " = H"(^f) , (33) 

is stored in the computer, the A K can be computed from equation (8). 

3.3 A Special Class of Weighting Functions 

The amount of computer memory required to store the set \H m \ 
will be small if h(t) is a function such that H m = for M < \m\ ^ N/2 
and M is a relatively small number. The H(co) corresponding to this class 
can be expressed by a particular form of equation (29) : 

U f, 

I Ufa) = Te~ iX sin X X) v . (34) 

n = -M -A — niT 

where 

X = uT/2 (35) 

and M <<C N/2. We have restricted our attention to the results that can 
be obtained with this class of weighting functions. 

Most of the well-known weighting functions, such as Hanning, 5 
Hamming, 5 and Taylor 10 ' 11 are in the class defined by equation (34). 
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The discrete Dolph-Chebychev and the Kaiser-Bessel 12 weighting 
functions, however, are not. 

The right side of equation (34) can be written over a common de- 
nominator to obtain the form 

H(») = Te~ iX S -^ - P(X) , (36) 

II (X 2 - nV) 

n = l 

where P(X) is, in general, a complex polynomial in X. We will restrict 
our attention to H(a>) with real C n and C = 1. In which case C_„ = 
C„ and, if D n = 2C n , we have 

h(t) = h r (t)[l + ED n cos(^)]- (37) 

Equation (34) becomes 



ff(co) = T*~" sin X 



l + f- M. 1. 

X £% X — n 7r J 



(38) 



In the next few sections wc will discuss three classes of weighting 
functions with the form of equations (37) and (38). They were chosen 
to provide two extreme cases of weighting and an intermediate example. 
Many other weighting functions in the class defined by equations (37) 
and (38) exist; the ones examined below provide sufficient data for our 
purposes. 

3.4 Class I Weighting Functions 

We first consider the class of weighting functions that provides the 
best possible reduction in \ H(u) \ for large w. Let this class be known 
as Class I. 

The only part of equation (36) that can be adjusted is the polynomial, 
P(X). Thus we must choose the coefficients, D n , to minimize | P(X) \ 
for large X. This is done by forcing P(X) to be a constant. The constant 
term in P(X), from equation (34), is 

P(0) = (-1) m t 2M (M\) 2 . (39) 

Hence, the desired class of weighting functions has the form, from 
equation (36), 



n* 2 - 



n ir 
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We denote the coefficients, D n , of this class of weighting functions 
as Di(M, n), making the dependence upon M explicit. From equations 
(36), (38), and (40) the Dj(M, n) are given by 

(-l)V Af (M!) 2 (X 2 -wV) 



D t (M,n) = lim 

rfli 2 - /<v 

Evaluation of the limit and some simplification gives 
2(-lY(Ml 



(41) 



D l (M,n) = 



■^B*^' (42) 



(M - n)\ (M +n)! 

We denote the weighting functions that use equation (41) as h M (t). 
Then from equation (37) 



hM = h T (t)\ 1 + £ Dj(M t n) cos 

L n = l 

This can easily be shown to be the same as 



2vnt \ 



(43) 




Fig. 2 — Spectra of the Class I weighting functions. 
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MO = MO 



4 Jtf (M!) s 
(2J1/) ! 



sin 1-=; 



(44) 



Thus the Class I weighting functions are described by equations (40) 
and (44). The so-called hanning weighting 5 is equivalent to hi(jt). Larsen 
and Singleton 13 used hi(t), h 2 (t), and others. 

Fig. 2 shows the shape of (l/T)e iuT/2 H M (a) for M up to 4. In Fig. 3 
we have plotted the normalized transmission of £Tir(w) (that is, 
20 Log 10 (H M (o))/T\). Several of the MO are shown in Fig. 4 and 
some values of Di(M, n) have been tabulated in Table I. 

3.5 Class II Weighting Functions (Taylor) 

Class I weighting functions provide the minimum high-order side- 
lobe amplitude in H(u) that is possible with a given value of M. We 
now turn to the class that gives the minimum main-lobe width, at the 
expense of higher side-lobe amplitude. 

The so-called continuous Dolph-Tchebycheff weighting functions 8 




-60 - 



-80 - 



-100 - 



-120 



\ \ WM=4 




Fig. 3 — Normalized loss of the Class I weighting functions. 
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Fig. 4 — Time response of the Class I weighting functions. 



provide the minimum main-lobe width in H(u>), consistent with a 
specified maximum side-lobe amplitude, but they are unrealizable 
functions. The Taylor 10 approximation to the Dolph-Tchebycheff func- 
tions provides almost the same main-lobe width and side lobes that have 
the specified maximum amplitude near the main lobe and then gradu- 
ally decrease as o increases. 

The Taylor functions have the form given by equations (37) and 
(38) with the D n dependent upon M and the maximum side-lobe ampli- 
tude, 1/R. We will denote the D n coefficients of Taylor weighting by 
D n (R, M, n), making the dependence upon R explicit. After adapting 
Taylor's equations to our situation, the D u 's are given by 



D u {R,M,n) = -■ 



fir.- 



in/ 



r + {K-\ 






R [' - (1)1 



(45) 



where 



and 



R = cosh (7rX) 



(M + l) 2 
X 2 + (M + \y 



(46) 



(47) 
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Table I — Values of D r (M, n) for M up to 4 





n 


M 


12 3 4 


1 
2 
3 
4 


-4/3 1/3 — — 
-3/2 3/5 -1/10 — 
-8/5 4/5 -8/35 1/35 



Solving equation (46) for X gives 

X = - In [R + VR 2 - 1]. (48) 

IT 

We will refer to the Taylor functions described by equations (45) 
through (48) as Class II weighting functions and denote them by 
k u (R, t) and K M (w). References 10 and 11 give a further discussion of 
Taylor functions. Taylor weighting functions have the property that, 
if M is too small, the D's given by equations (45) will define an H(w) 
whose first few side lobes have the amplitude given by equation (46), 
but some of the higher-order side lobes will have much higher ampli- 
tudes. Thus, for each value of desired side-lobe level, l/R, there is a 
minimum value of M that will give good side-lobe suppression. 

Some minimum values of M that give good side-lobe control are 
listed in Table II. 

Figures 5 and 6 show the shapes of a Class II weighting function with 
M = 7 and R = 10 3 . This particular weighting function will be ex- 
amined below when simulation results are compared. It will be shown 
there that this weighting function is useful when the received tone 
frequencies are very closely spaced. 

3.6 Class III Weighting Functions 

The third class of weighting functions has been chosen to have, to a 

Table II — Minimum Values of M for Good Side-Lobe Control 



20 log,„ R (db) 


M 


36 


3 


42 


4 


48 


5 


54 


6 


60 


7 


66 


9 
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-20 - 



-40 - 




-60 - 



Fig. 5 — Normalized loss of the Class II weighting functions. 

large extent, the desirable properties of both Class I and Class II 
weighting functions. That is, Class III weighting provides better reso- 
lution than Class I weighting for tones with a "small" frequency 
separation. Moreover, they also provide better resolution than Class II 
weighting of tones with a "large" frequency separation. 

We will identify the Class III weighting functions by g M (t) and 
GmW), where g M (t) <-> G M (u>). The Z>„ coefficients for this class will be 
denoted by D U \{M, n). The first member of the class is chosen as dis- 




Fig. 6 — Time response of the Class II weighting functions; M = 7, R = 1000. 
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cussed below and the other members are obtained by operations on the 
first. 

In order to have the high-order side lobes of Gm(u>) fall off at least 
as 1/co 2 we must have g M (0-{-) = g M (T—) = 0. In terms of the D„ 
values this means that 

1 + Z-Din(M,n) = 0. (49) 

n=l 

With this restriction, of course, gi(t) is the same as the Class I weighting 
function, h x (f). Thus the distinguishing properties of Class III weighting 
are determined by g 2 (t). 

We chose the coefficients of g 2 (t) so that the loss of G 2 (o)) reached 60 
dB with as small a value of co as possible with the side lobes of G 2 (u>) 
never exceeding —60 dB, subject to equation (49). The D n values for 
this condition are: 

D m (2, 1) = -1.19685, D IU (2, 2) = 0.19685. 

The g 2 {t) thus defined is almost the same as Blackman's 8 proposed 
function, Q 4 (/)- 

The rest of the members of the Class III functions are defined in a 
manner similar to that used by Helms 6 for the synthesis of digital 
filters. We define 

(f . h M - 2 (t)g 2 (t) 

9uW 1 + * Z) in (2, 1) Dy{M - 2, 1) + | D,„(2, 2) D 1 (M - 2, 2) ' 

M > 2. (50) 

The normalization in equation (50) puts g M (t) in the form of equation 

< 37) - 
The Class III functions just defined have high-order side lobes, in 

G M (oi), that decrease as u~ M . This contrasts with oT (M+1> for Class I 

weighting and with co -1 for Class II weighting. Thus, Class III weighting 

functions provide slightly narrower main-lobe width than Class I at the 

expense of slightly higher side lobes. 

Some values of Dm(M, n) are tabulated in Table III. 

In Fig. 7 we have plotted the normalized spectra (that is, 
(l/T)e iaT/2 G M (a>)) of some Class III weighting functions. Fig. 8 illus- 
trates the normalized loss provided by (? w (w) for values of M up to 4. 
It is interesting to note, from Fig. 7, that G 2 (o>) reaches —50 dB before 
any of the others, just as H 2 (co) did in Fig. 3. In Fig. 9 we have plotted 
g M (t) for values of M up to 4. 
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Table III — Values of D iu (M, n) for M up to 4 





n 


M 


1 


2 


3 


4 


2 
3 
4 


-1.19685 
-1.43596 
-1.566272 


. 19685 

.497537 
.725448 


-.0615762 
-.180645 


.0179211 



rv. response to a cosine wave 



We are interested in measuring the frequencies and levels of signals 

that comprise several sine waves. In view of this and the linearity of the 
DFT it is convenient to examine the properties of the DFT of samples 
of cos (o) t + 0). 



4.1 Basic Formulas 
Let 



s(t) = cos (fti t + 0) 

f(t) = h T (t) cos (a> t + 0). 



(51) 
(52) 




Fig. 7 — Spectra of the Class III weighting functions. 
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Fig. 8 — Normalized loss of the Class III weighting functions. 
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Fig. 9 — Time response of the Class III weighting functions. 
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and 

H T (u) = f h T (t)e~ iat dt (53) 

or, in the usual notation, 

H T (co) <-> h T (t). (54) 

Then 

F'(») - §e"ff£(cu - co ) + !e-"H*( W + co ) (55) 

and from equation (19) 
A(«) = ie"H' T (a - co ) + |e-'''H*(w + w ) 

+ §[cos 6 - cos (wo? 1 + e)e~ iuT ]. (56) 
With A r (i) as denned by equation (22) 

ffr(B )- re -'"- S -S^. (57) 

It is possible to put equation (57) in equation (56) and evaluate the 
indicated summations. The same answer can, however, be found more 
easily from equation (18). With equation (52) in equation (18) we have 

a(«) = z cos (^f + 9 y*»r,* (58) 

or 

Aico) = y Z e- i{a -°' )nT/N + y" E e - llm * mM '". (59) 

n = n = 

Both sums in equation (59) are finite geometric series. Thus after re- 
arrangement we obtain 

A ( U ) = y v"*-"- sJ ^ + i.-'V"*-"" sJ ^ (60) 

v ' l sinz smy 

where 

2 = ^i^ < 61 > 



and 



» - Hf^- (62) 
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The formulation given by equation (56) is important from a con- 
ceptual point of view while that of equation (60) is useful for numerical 
evaluations. 

The evaluation of equation (60) when either z or y is zero requires an 
examination of limits, which can be done by inspection. 

4.2 General A(u>) 

The use of one of the weighting functions denned by equations (34), 
(37), and (38) gives rise to an A(<a) different from the one calculated in 
equation (60) . The more generalized form of A (co) is given by 

U Hz-nir/N) 



A ( W ) = y e e~ iN ' sin Nz £ C» 



-tt . I nir 

sm ri 

M Hv-nr/N) 

+ K'V'^sinM/ £ C-^ r- (63) 

n — u . I nir 



sin \y - jj 

All of the simulations, to be discussed below, used this A(o>), in 
equation (21), to compute .A* values. 

4.3 Approximations 

We will make use of several approximations in the next section. The 
important ones are established here. In this section h T (t) is assumed to 
be the weighting function and co is in the range < co < w,/2. 

Consider equation (19). If | F(w — he.) | is small for I ^ 0; then 

A(a>) ^ j- F(«) + M/(0+) - KT-)e~ iaT ] (64) 

where 

t, = T/N. (65) 

Thus, from equation (21), 

A K ^ |- f(^) + i[/(0+) - f(T-)]. (66) 

Next consider equations (56) and (57). | H T (o>) | is "large" only 
near w = 0. Thus we obtain another approximation, used when s(t) = 
cos (o} t + 0), 

F(w) « y $ H T (u - uo). (67) 
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From equation (57) we then obtain 



1 N ,-,,-,-| (B - M .)r/2i sin [(co - u a )T/2] 

l F ^=^ 6e [(co - v )T/2] 



(68) 



Thus, 



Ajc | = 2 



sin VW ~ ~2 



or, because, 



Ku.T 

2N 



Kw,T _ Wo.T 

2iV 2 



= Ktt, 



(69) 



A l~5 



sin(&-^ 



7<L7T - 



2 



(70) 



(71) 



From equation (71) we see that, apart from the error in the approxi- 
mations, one should be able to accurately estimate the frequency and 
magnitude of a cosine wave from the A K values for K near o T/2t. 

The main lobe of a sin X/X curve reaches zero at X = ±t. Thus the 
main lobe of the curve, of which the A K are points, reaches zero at 
o> = a>o ± 2ir/T or, since <a,/N = 2tt/T, at 

co = cu rfc k),/N '. 

The main lobe is just wide enough to contain two A K values (estimators), 
except if co equals some multiple of u,/N. It will be shown later that 
two A K values will be enough to estimate the parameters of a cosine 
wave. 



V. FREQUENCY AND LEVEL ESTIMATION 

In the preceding sections we have developed methods of computing 
the DFT coefficients of a known input (for simulations) and have dis- 
cussed three classes of weighting functions to improve the effective 
selectivity of the DFT process. Our final task, which is undertaken in 
this section, is to determine accurate ways of processing the DFT 
coefficients to extract the frequencies and magnitudes present in the 
sampled signal. 

The methods are to be useful when the real-time demands upon the 
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computer are important. Thus all of the methods were chosen to have 
simple formulas which the computer can be programmed to evaluate 
when it is making a measurement. 

The equations we have examined fall into two classes: those that use 
only two estimators to calculate the frequency and level of a cosine 
wave and those that use many. The following paragraphs derive the 
most promising of these equations. The next section will present the 
accuracies that the various equations can achieve. We start with a 
formula that makes use of many estimators. 

5.1 Method 1 

The derivation of this method is somewhat involved, so we first 
explain the motivation behind it. 

Suppose one has an f(t) that is known to be given by 

/(0 = B cos o> Q t, (72) 

and one wants to determine w and B from operations on j(t). One way 
to determine B is from 

limi f T f(t)dt= limi f B*/2dt 

= B-/2 (73) 

thus, 

B 2 = 2 1im^ [ T f(t)dt. (74) 

The derivative of jit) is 

f(t) = -Bu smu t (75) 

and 

lim^ f ' [j'(t)fdt =BW /2 (76) 

Thus w can be determined from 

lim^ f [j'{t)] 2 dt 

"o = —77 ( 77 ) 

lim ± \ [jit)f dt 

r—oo J ■'o 

The above result is the motivation for the following derivation. 
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Assume 

s(t) = B cos M + 6) (78) 

and 

/(0 = h(t)B cos (co « + 0) (79) 

where h(t) is one of the weighting functions, given by 

HD - [l + E A. cos (^)]m0- (80) 

Let the estimators, A K , be given by equation (31). 
Define 

• T 



Po-lim^f [/(fl] 8 dl (81) 



and 



P, - lim^ [ T [fiQVdt. (82) 

T-^co •* •'0 

Expansion of equation (79) gives 
/(*) = Bh T (t)\ cos (o3 t + B) 

+ \ E A.[cos ("o* - w« i + ») + cos (w * + na> a * + *)] \ (83) 

^ n-l J 

where 

co a = 2tt/^. (84) 

From equations (81) and (83) 



Po = f{i + |E4 



(85) 



Calculation of Pi from equation (83) yields, assuming nu a j£ o> , 

d2 ( m rfi 2 Af ,»» 2 n 2 , . 2 1 

P. - f {«! + Z ^f + E HH' t 88 ) 

Combining equations (85) and (86) gives 



2 Pi 

"0 = p" 

* 



2+ E^n 



(87) 
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The next step is to compute P,/P„ from the DFT of j(t). Using the 
approximation 



J-V k--N/2 



(88) 



which can be rearranged to give 

N/2 
11 1 .1 » 

N 



j(t) w ^ <U + 2 £ [Re (A K ) cos Kco * - Im (A K ) sin £»„*] 



From equation (89) 



and 



Thus, 



1 Nn 



1 iV/2 

*" IC-i 






(89) 



(90) 



(91) 



p <2 AV2 

^ K-l 



(92) 



We use equation (92) in equation (87) to obtain the final result. 
Denote the estimated co by co . Then 



a2 2J 

OJ„ = CO 



E # 2 1 A, 



-V/2 

E 

K=l 



- u. 



Ul + Z I A, 



(93) 



where 



E»'tf 



c/m = — 



(94) 



2+ £# 



By Z>„ we mean, of course, D^M, n), D U (R, M, n) or D ul (M, n). Thus 
Method 1 is applicable to all three classes of weighting functions. Some 
values for U u are tabulated in Table IV. 
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Table IV — Values of U m and V m 





Clasa I Weighting 




M 


U M 


v M 



1 
2 
3 

4 



.333333 
.571429 

.818182 1 
1.06667 1 

Class II Weighting 

(M = 7, R = 1000) 


.5 

.75 

.972222 

.155 

.31327 


U-, 


= .357071, V 7 = .769710 




Class III Weighting 




M 


U M 


V M 


2 
3 

4 


.45732 

.715523 1 
.968949 1 


.8678 

.07833 

.25033 



The way to use equation (93), when more than one tone is present 
in s(t), is to use only the estimators, A K , for K tt co /co a to calculate 
each a> • The simulations below will show that this technique gives 
accurate results. 

5.1.1 Estimation of Level 

From equations (85) and (90), we obtain the way to estimate B. 

N/2 

, A\ + 2 £ | A K | 2 

B ~ W v„ 

where 



(95) 



Vm = h + \ E & • (96) 

n = l 

Some values of 7 A / are tabulated in Table IV. Observe that if only 
the basic weighting, h T (t), is used, then U M and V M become zero and 
|, respectively. 

5.2 Method 2 

The preceding formulas for estimating the frequency and level of a 
cosine wave from its DFT use more than two estimators. The next few 
paragraphs establish formulas that require only two estimators. The 
formulas apply only to the Class I weighting functions, h M (t). 
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We start by recalling the approximations given by equations (64) 
and (67). With H M {u>) substituted for // r (co) in equation (67) we have 



A(u)tt^-e i6 H M (a -«o). 



From equation (40) 



//, 



(2X\ _ _,, v sjnX (-DV^Mif 



\T 



= Te 



X 



II (* 2 - nV) 



(97) 



(98) 



where X = coT/2. If s(t) = B cos (o) t + 0) then, from equations (97) and 
(98), 



A(a 



„ , , lN2 13JV sin7rfl 

9 * 



:-i)- v (M!) 2 ^ 



t II (»+n) 



where 



i' = (co — u)/a} a . 



(99) 



(100) 



Suppose the largest* estimator is A , and its largest immediate neighbor 
A,. . Of course m — I = ±1. Define 



Let 



and 



Define 



a = m — I = ±1. 

a, = | A t | 

a 2 = | -4 m | = | Ai +a !■ 

u = ^ _ j. _i/2 ^ w ^ 1/2. 



Then from equation (99) 

o, w(-i)-(af!r 



7T n (« + ») 



(101) 
(102) 
(103) 

(104) 
(105) 



t By largest we mean | A\ | > | ^4ir | for K ^ /. 
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and 

a, w (-!)"(«!)* ^ r' (M ~ a) (106) 

T IJ (w + n — a) 

n=-AT 

since a = ±1, equation (106) is the same as 

a, « -(-.)'(«' f — ^ (107) 

x n (« + ») 

n= — Af— a 

Division of equation (105) by equation (107) gives 

<h ~ u — 2JM. + 1) 

a 2 ~ w + «M 

Define 

a 2 (M + 1) - aiM 

M x = : 

a, + a 2 
then from equation (108) an estimate of u is 

u = aWi . (110) 

Hence, the estimate of co is given by 

do = *> a (l + w). (Ill) 

From equation (105) the estimate of B is 

M 



(108) 



(109) 



a 1 2x(-l) Af II (A - n) 

5 = N{M^l\nra » (U2) 

where li is defined by equation (110). 

Another version of equation (112), better for machine computation, is 

*- J^f 5? ftl -(-)'• ( 113 ) 

AT sm (ttw) „=i W 

Method 2 with M = is essentially the same as was derived by 
Penhune and Martin 14 to solve a radar problem. 

5.3 Methods 

The simplicity of the estimation equations of Method 2 led us to 
extend this method to include any class of weighting functions described, 
in general form, by equations (37) and (38). We will refer to this more 
general method as Method 3. 
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From equations (109) and (110) we see that if a Class I weighting 
function is used, one way to obtain an estimate of u is to use the function 

Ea 2 + a x 

in equation (110). There are three degrees of freedom in the bilinear 
form and we chose to express them in the manner shown in equation 
(114). 

Method 3 is simply the application of equation (114) to other classes 
of weighting functions. We obtained values for the coefficients in equa- 
tion (114), for several weighting functions by: 

(i) Computing a, and a 2 , using equations (21), (63), (102), and 
(103), with 9 = and many values for co near w,/4. 

(ii) Computing the corresponding values of u from equation (104). 
(iii) Choosing values for C, D, and E such that equation (114) gave 
a good fit to the data computed in the first two steps. It turns out that 
the curve described by equation (114) is satisfactory if it fits the com- 
puted data exactly at u = 0, j, and J. 

In this manner we obtained the following coefficient values: 

Class II weighting, M = 7, R = 1000; 
C = 1.96339, D = 1.01643, E = 0.893534. 

Class III weighting, M = 2; C = 2.56919, 
D = 1.5374, E = 1.06345. For M = 3; 
C = 3.6020, D = 2.5862, E = 1.0317. 

Using an approximation similar to equation (71), but extended to 
include weighting functions described by equation (38), we obtain an 
estimate for B, 

6 27rai 

(115) 



w^+f;^] 



Method 3 can be used with weighting functions specified only in 
terms of H* (a>) as well as those given in terms of #(aj). 



VI. COMPARISON OF ACCURACIES 



In the preceding paragraphs we have derived several formulas that 
produce estimates of w and B from A K values (estimators). We now 
turn to a comparison of these formulas on the basis of accuracy. The 
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accuracies we will compare do not include any possible computation or 
nonlinearity errors or other accuracy limitations that may be present 
in a DFT analysis system. Our accuracy comparisons include only the 
effects of leakage. 

The estimators used in the simulations were generated by using the 
function A(<a), described by equation (63), in equation (21). This is 
equivalent to applying the weighting by multiplication in the time do- 
main or by convolution in the frequency domain. The use of equation 
(63) in simulations greatly reduces computation time. All of the simula- 
tions used N = 512 and /. = N/T = 7040 Hz. 

All of the estimation methods presented above use approximations. 
In this section we shall demonstrate just how good the approximations 
are. 

Consider the case where a tone of frequency / , angle O , and ampli- 
tude B is being measured while another tone, at frequency /i , angle 
0i , and level B t , is also being received. The presence of /i will affect the 
accuracy of any estimate one makes of / or B (due to leakage). The 
size of the errors in the estimates of / and B will depend upon which 
formula (method) is used and upon the values of / , B n , 0„ , ji , B v , 
and 0i . The combination of parameters that causes one method to 
give the worst estimates will, in general, not be the combination that 
causes another method to be at its worst. Thus it is difficult to compare 
methods. 

We have compared the three methods on the basis of the worst 
estimates each will make when O and / are confined to a specified 
range of values (for example, 990 ^ /„ ^ 1003.75 Hz and ^ O ^ 360 
degrees) and B , B x , j l , and 0i are fixed (for example, B = B t - I 
and 0, = degrees). 

Notice that if /i is equal to some multiple of 1/T then its A K will be 
very small except for K near Tfi . Thus such an ] x cannot cause much 
error in any of the three methods of estimation of / , which use A K 
values. For this reason we have fixed ji at a value that is an odd multiple 
of |T. 

Tables V and VI illustrate how inaccurate frequency and level esti- 

Table V — Poorest Estimates with Interference Separation in 
the Range 55 ± 6.88 Hz, No Special Weighting 

Method Frequency Error, Hz Magnitude Error, dB 

It 6.82 .596 

2 7.37 .580 

t Method 1 using only six estimators, those from I — 2a to / + 3a. 



DISCRETE FOURIER TRANSFORM 225 

Table VI — Poorest Estimates with Interference Separation in 
the Range 178.96 ± 6.88 Hz, No Special Weighting 

Method Frequency Error, Hz Magnitude Error, dB 

It 3.08 .0628 

2 3.75 .165 

t Method 1 using only six estimators. 

mates are when no special weighting (M = 0) is used. The interfering 
tones corresponding to the simulations described in Tables V and VI 
were located at 1051.88 and 1175.63 Hz respectively. The frequency 
and magnitude error entries in these, and subsequent, tables give 
absolute values only. 

The simulation results presented in Tables VII and VIII indicate that 
substantial improvement in the accuracy of frequency and magnitude 
estimates can be achieved when weighting is used. The data in Table 
VII shows that when the two tones are separated by a "small" frequency 
difference accurate frequency and level estimates can be obtained by 
using Method 3 with Class II or Class III weighting. Table VIII in- 
dicates that as the frequency separation increases Method 2 with Class I 
weighting is better. The accuracy of estimates made on closely spaced 
tones can, of course, always be improved by increasing N and T while 
keeping the ratio N/T constant. 

It is interesting to note from Figs. 3 and 8 that, for a given value of 

Table VII — Poorest Estimates with Interference Separation in 
the Range 55 ± 6.88 Hz 



Class I Weighting 










Method 


M 


Frequency Error, Hz 


Magnitude Error, 


dB 


1 


l 


.89 


.11 




1 


2 


3.32 


.47 




1 


3 


6.16 


.96 




2 


1 


.513 


.12 




2 


2 


.104 


.11 




2 


3 


.409 


.5 




Class II Weighting, 


R = 


1000 






Method 


M 


Frequency Error, Hz 


Magnitude Error, 


dB 


1 


7 


1.14 


.120 




3 


7 


.0651 


9.35E-3 




Class III Weighting 










Method 


M 


Frequency Error, Hz 


Magnitude Error, 


dB 


1 


2 


2.11 


.225 




1 


3 


5.12 


.821 




3 


2 


.034 


.026 




3 


3 


.149 


.0655 
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Table VIII — Poorest Estimates with Intererence Separation in 
the Range 178.96 ± 6.88 Hz 



Class I Weighting 










Method 


M 


Frequency Error, 


Hz 


Magnitude Error, dB 


1 


1 


1.67E-3 




2.79E-4 


1 


2 


3.98E-4 




5.25E-5 


1 


3 


4.55E-2 




4.70E-3 


2 


1 


5.73E-3 




1.42E-3 


2 


2 


1.47E-4 




2.19E-5 


2 


3 


2 . 13E-5 




1.55E-6 



Class II Weighting, R = 1000 

Method M Frequency Error, Hz Magnitude Error, dB 

1 7 5.56E-3 2.05E-5 

3 7 8.09E-2 1.07E-2 

Class III Weighting 

Method M Frequency Error, Hz Magnitude Error, dB 

1 2 9.10E-5 1.19E-5 

1 3 1.85E-2 1.90E-3 

3 2 3.35E-3 3.24E-3 

3 3 5.94E-4 9.62E-6 



M, there is not a great deal of difference between the weighting con- 
tributed by H M {u) or G M (a)). However, from Table VII it is obvious that 
the use of G a (<a), when the tones are close together, will yield much more 
accurate estimates than Class I weighting. 

In equation (101) the "pointer", a, was defined. The value of a is 
used by the system to determine whether the frequency being measured 
is above or below the frequency of the largest estimator, A t . Our 
simulation studies showed that under certain circumstances a, as 
calculated by equation (101) will point in the wrong direction. In 
general this happens when the contributions to A,_i and A i + l due to 
the interference is equal to or greater than the difference in the con- 
tributions to the same estimators due to the tone being measured. Thus 
if | Ai-i | « | Ai+i |, a small difference in their magnitudes can change 
a. In our simulations this effect only caused trouble when / and /i were 
separated by less than half the width of the main lobe of the weighting 
function, H(co). 

Since we have fixed B = Bi = 1 in the simulations we have ignored 
the adverse effects of "large" level differences on the accuracies of the 
various methods. Leakage from an interfering tone with a high level, 
relative to the tone of interest, would certainly tend to reduce the 
accuracy provided by any of the three methods, no matter which weight- 
ing is used. 
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VII. CONCLUSIONS 



The discrete-tone measurement system we have been discussing is 
particularly well suited to systems that involve computer-controlled test- 
ing or measurement, provided the real time needed for the computations 
is available. Two advantages are: 

(i) The only interface hardware is the A-D converter (with its 
lowpass filter). 

(ii) The system is capable of measuring many received frequencies 
and levels during the same computation time. 

For a given number of samples taken at a given sampling rate, the 
accuracy of the system can be significantly improved through the use 
of some type of weighting function. We have examined three classes of 
such smoothing functions and have developed formulas which permit 
the extraction of received signal frequencies and levels from the DFT 
coefficients. The results of system simulations, presented in Tables V 
through VIII, show that the inherent accuracy of the described system 
can, through the proper use of weighting functions and estimation 
methods, be made satisfactory for many applications. 

With Method 2 considered to be a special case of Method 3, the tables 
show that the best estimation method, for all of the weighting functions 
examined, is Method 3. 

The tables also show that there is no "best" weighting function. The 
weighting to be used for any particular application should be selected 
only after a consideration of the expected tone frequencies, the relative 
levels, the measurement accuracy desired, and the desired value for N. 
The sampling frequency, N/T, should be more than twice the highest 
frequency to be measured. 

It is interesting to observe that the Taylor (Class II) weighting func- 
tion used in the simulation is, for the situations simulated, not signifi- 
cantly better than the Class III weighting, G 2 (t), which is essentially 
that proposed by Blackman. 5 There may be other situations, however, 
when the nearly optimum main-lobe width of the Taylor functions is 
useful. 

If the system could tolerate the relatively large amount of computer 
memory required, then the discrete Dolph-Chebychev functions de- 
scribed by equation (27) could provide some advantages. 
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